Supplementary Report for the article Functional diversity can facilitate the collapse of an undesirable ecosystem state

1 Introduction

This document is the Supplementary Report for the report Functional diversity can facilitate the collapse of an undesirable ecosystem state.

2 Creating among strain diversity

Here we provide additional description of the method used to create strain variation within functional groups. The following should be read only after the relevant text in the main report. We generated variation in maximum growth rate \(G_{max}\) and tolerance to reduced sulfur \(h_{SR}\) or oxygen \(h_{O}\) by varying the range of trait values around a reference trait value. To manipulate the amount of variation in a given trait and functional group, we introduced the meta-parameter \(Div\) subscripted by the trait and the functional group. Take for example the cyanobacteria functional group, a reference maximum growth rate parameter \(G_{max,CB}\) and a reference tolerance parameter \(h_{SR,CB}\). The amount of among strain variation in the cyanobacteria is controlled by the two meta-parameters \(Div_{G_{max,CB}}\) and \(Div_{h_{SR,CB}}\). To create a trade-off, the two meta-parameters were always of different sign. The among strain variation was calculated such that the growth rate of strain \(i = 1\) was a factor \(1/({2Div_{G_{max,CB}}})\) of the reference growth rate, and the growth rate of the \(i = n\) strain was a factor \(2Div_{G_{max,CB}}\) of the reference growth rate. The growth rate of the \(i = {2,...n-1}\) other strains was the reference growth rate multiplied by a factor that was equally distributed between the factor of the \(i = 1\) and \(i = n\) strain. This implementation of the variation and trade-off has the assumption of a linear trade-off of log-log transformed among-strain variation. Hence, for example, with the reference growth rate \(G_{max,CB} = 0.05\), among strain variation of \(Div_{G_{max,CB}} = 1\), and nine strains (\(n = 9\)), the growth rates of the nine strains are 0.025, 0.03, 0.035, 0.042, 0.05, 0.059, 0.071, 0.084, 0.1. Put another way, the growth rate parameter of the \(i\)-th of the \(n\) strains in a functional group was calculated as \(G_{max,i} = D_i * G_{max}\) where \(D_i\) is the \(i\)-th element of the set \(D = {2f(x)Div_{G_{max}} |x ∈ N, x ≤ n}\), where \(f(x) = (x − (n + 1)/2)/((n − 1)/2)\).

This procedure ensures that the range of trait values is independent of the number of strains. When there were only two strains, the two strains took the minimum and maximum traits values, as defined by the \(Div\) variable for that trait. When there were three strains, one had the minimum trait value, one the maximum, and one the mean. With nine strains, one had the minimum trait value, one the maximum, one the mean, and three were equally spaced (on log scale) between the minimum and the mean, and the remaining three between the maximum and the mean.

We show three examples of the created variation among strains within functional groups, one for nine strains (Figure 2.1), one for three (Figure 2.2), and one for two (Figure 2.3). In each of these three examples, three graphs show the diversity among strains in the maximum growth rate trait and the tolerance trait, plotted against each other and therefore also showing the trade-off between the two traits.

Notice that in all three cases, the range of trait variation represented by the strains is the same, even though the number of strains differs. This approach allowed us to independently manipulate trait diversity and number of strains. Also note that maximum growth rate of strains is evenly distributed on a log scale, even though it may appear at first site that they are distributed evenly on a linear scale (e.g. Figure 2.2).

All three examples have the following amounts of trait variation: \(Div_{G_{max,CB}}\) = 0.032, \(Div_{h_{SR,CB}}\) = -0.16, \(Div_{G_{max,SB}}\) = 0.095, \(Div_{h_{O,SB}}\) = -1.938, \(Div_{G_{max,PB}}\) = 0.095, \(Div_{h_{O,PB}}\) = -1.938. These are also the maximum amounts of trait variation explored in the simulations.

Figure 2.1: Trait variation with nine strains in each of the functional groups. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

Figure 2.2: Trait variation with three strains in each of the functional groups. Note that the three displayed strains are the 1st, 5th, and 9th strains from the nine strain example shown in Figure 2.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

Figure 2.3: Trait variation with two strains in each of the functional groups. Note that the two displayed strains are the 1st and 9th strains from the nine strain example shown in Figure 2.1. In all panels colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains. Tolerance is given as the concentration of the inhibiting substrate where the growth rate of the strain is reduced by 50%.

A consequence of this variation and of the trade-offs is shown in Figure 2.4, which shows the relationship between realised growth rate and the concentration in the environment of the respective inhibiting substance, e.g. reduced sulfur for the cyanobacteria. In these graphs only the least and most tolerant strains are shown. The lines cross due to the trade-off. At high concentrations of the inhibiting substance the more tolerant strain has higher realised growth rate, while at lower concentrations it is the less tolerant strain that has the higher realised growth rate.

Figure 2.4: Response of realised growth rate to variation in concentration of the inhibiting substrate for two strains, namely the most and the least tolerant strain. As always colour indicates the functional group, darker shades show more tolerant-lower growth rate strains, lighter shades show less tolerant-higher growth rate strains.

3 Method for finding stable states

3.1 General method

We now discuss two approaches for finding how the stable state of the system varies along the gradient of oxygen diffusivity when there exists the possibility of multiple stable states. In the main article we describe and report results of the method that we term the Temporal method.

Here we also describe and give results of what we term the Replication method. In this, we made separate simulations for each of many values of oxygen diffusivity, ran the simulations for long enough to ensure transients had passed, and recorded the state. We ran multiple simulations per value of oxygen diffusivity, and started them with different initial conditions (either low or high total cyanobacterial abundance), so as to check for presence of alternate stable states. This method was used by Bush et al 2017. We term it the replication method since it is achieved by having many independent replicates of the ecosystem, each with a different oxygen diffusivity, and also a separate replicate for each of the different chosen initial conditions.

As is shown in the figures below, the temporal method results in a larger region of bistability than the replication method, but only during increasing oxygen diffusivities. The difference in results occurs because, as discussed below, the initial conditions differ between the two approaches.

3.1.1 Replication method, no diversity

Figure 3.1: Stable states found with the replication method. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

3.1.2 Temporal method, no diversity

Figure 3.2: Stable states found with the temporal method. All other aspects are the same as in Figure 3.1.

3.1.3 Both methods overlaid, no diversity

Figure 3.3: Stable states of the two methods overlaid, the temporal method with a dashed line, and the replication method with a solid line.

3.2 Effect of diversity, replication method

3.2.1 Overall effects of diversity

The replication method shows the same effects of functional diversity as those revealed by the temporal method (Figure 3.4). Namely

• Introducing functional diversity in only the cyanobacteria increases the resilience of the oxic state but has no effect on the resilience of the anoxic state.
• Introducing functional diversity in only the sulfate-reducing bacteria increases the resilience of the anoxic state and reduces the resilience of the oxic state.
• Introducing functional diversity in only the phototrophic sulfur bacteria reduces the resilience of the anoxic state and has no effect on the resilience of the oxic state.
• Introducing functional diversity in more than one functional group reduces, reverses, or does not change the effect of diversity in individual functional groups. For example, introducing functional diversity in all three functional groups reduces the resilience of the oxic state when trait variation is low but increases the resilience of the oxic state when trait variation is high.

Figure 3.4: Replication method results. Left column of panels: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. Right column of panels: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience.

3.2.2 An example of stable states

Figure 3.5 shows an example of stable states along the oxygen diffusivity gradient produced by the replication method for finding stable states. The patterns are similar to those produced by the temporal method for finding the stable states (see Figure 4 in the main report), though they differ in one respect: close to the transition back to the state a group dominates, some of the more tolerant strains occur in the stable states found by the temporal method (Figure 4 in the main report), while they do not occur in stable states found by the replication method (Figure 3.5 here). We understand that the appearance of these strains is transient, and the difference in their appearance results from the differences in starting conditions of the two methods. Specifically, in the temporal method, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the replication method the starting conditions are identical at each level of oxygen diffusivity. These differences in starting conditions apparently result in longer persistence of transients in the temporal method than the replication method.

An interesting feature of these results is that although the position of the oxic to anoxic tipping point is affected by diversity in the sulfate-reducing bacteria (Figure 3.4 and Figure 3 in the main report) there is no evidence of more tolerant strains of sulfate-reducing bacteria being present in the final system state (Figure 3.5). This phenomenon is explained in Section 5.2.

Figure 3.5: Stable states found with the replication method with intermediate levels of diversity in all functional groups. Figure features are otherwise the same as those in Figure 3.1).

3.3 Effect of abundance change event

As mentioned in the main report, to avoid very low strain densities, we added 1 unit of density to every biological state variable every 1,000 hours. We here show that a slightly different method produced the same results. Specifically, we tested the effect of making an abundance floor of 1 unit that was implemented every 1,000 hours. I.e. every 1,000 hours, any abundances that were lower than 1 unit were set to 1 unit. As shown in Figure 3.6, the results with this floor abundance are similar as when 1 was added every 1,000 hours. Note that a coarser gradient of trait variation was used, hence the saw-tooth patterns. A further difference is the large decrease in resilience of the anoxic to oxic transition above moderate levels of diversity in the SB and CB-SB treatments. This pattern can also occur with addition of 1 abundance unit (e.g., see Figure 4.3), and is discussed in the section that concerns this below (section SB diversity effects).

Figure 3.6: Results with an abundance floor. Left column of panels: the effect of trait diversity (y-axis) on the positions of the tipping points (anoxic to oxic in blue, oxic to anoxic in red), which determine the extent of the region of bistability (green shaded region) in the oxygen diffusivity gradient. The left-most region is anoxic only, the middle is oxic or anoxic depending on historical conditions, and the right-most is oxic only. The thick grey vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations. Right column of panels: The effect size of trait variation on resilience, with positive values on the y-axis indicating that diversity has increased resilience, i.e. delayed the shift to the other state.

3.4 Further thoughts

Here we give some further thoughts regarding the comparison of the temporal and replication method for stable state finding.

A particularly notable difference between the two approaches is in the initial condition. In the temporal approach, the initial condition at a new level of oxygen diffusivity is the stable state at the previous level of oxygen diffusivity, while in the approach of Bush et al 2017 the initial conditions are fixed and arbitrary (though this approach is still useful for illustrating the bistability of the system). This difference in initial conditions increases the region of bistability in the temporal approach because as a tipping point approaches, the functional group that will dominate after the regime shift starts from much lower abundances than in the approach by Bush et al 2017 and thus requires stronger amelioration of the environment to take over. Nevertheless, our main conclusions about the effects of functional trait diversity are the same using either approach.

The temporal approach is particularly useful when investigating the effects of strain richness. In preliminary simulations using the Bush et al 2017 approach, we found the counterintuitive result that higher strain richness led to lower resilience. We then recalled that we chose to divide the initial total abundance of the cyanobacterial (CB) functional group equally among the strains (termed a substitutive design in biodiversity manipulation experiments). This meant that the most tolerant CB strain in a three-strain simulation had higher initial abundance than the same most tolerant CB strain in a nine-strain simulation. The greater abundance of the most CB tolerant strain made it possible for a three-strain system to develop to be oxic, whereas the nine-strain system developed to be anoxic, despite identical range of trait variation. We could have chosen to do otherwise, for example start with the least tolerant strains of cyanobacteria as more abundant when the system started with high cyanobacterial abundance, and the most tolerant strain as more abundant when the system starts with low cyanobacterial abundance. Using the temporal approach removed the need to specify initial abundances (except the very first one), and therefore results were not so determined by the assumptions used to determine the initial conditions.

4 Effect of simulation duration

4.1 Diversity effects

Recall that the temporal method for finding stable states involves running the simulation for a certain amount of time, termed the “wait time”, at a constant level of oxygen diffusivity. At the end of that wait time the state of the system is recorded, and then the oxygen diffusivity is slightly increased or decreased, and the system then simulated again for the duration of the wait time.

The following figures 4.1-4.5 show the effects of diversity with long to shorter wait times.

Figure 4.1: Effects of diversity with wait time of 1e6 hours (same figure as in main report). (All else as in (Figure 3.4).

Figure 4.2: Effects of diversity with wait time of 1e5 hours. (All else as in (Figure 3.4).

Figure 4.3: Effects of diversity with wait time of 1e4 hours. (All else as in (Figure 3.4).

Figure 4.4: Effects of diversity with wait time of 1e3 hours. (All else as in Figure 3.4).

Figure 4.5: Effects of diversity with wait time of 1e2 hours. (All else as in Figure 3.4).

4.2 Example of transient dynamics at 1e6 hours

The results below show the stable states as a function of oxygen diffusivity for a system with an intermediate level of diversity in all functional groups, with either 1’000’000 hours at each oxygen diffusivity level (Figure 4.6), or 2’000’000 hours (Figure 4.7). This shows that 1’000’000 hours duration of the simulations are sufficient to reach stable state, except very close to state transitions on the trajectory of recovery (i.e.: cyanobacteria: trajectory of increasing oxygen diffusivity; sulfate-reducing bacteria: trajectory of decreasing oxygen diffusivit). Very close to these transitions, even at 2’000’000 hours, the dynamics are still transient, with high tolerance strains at high density, even though they will eventually be outcompeted by low tolerance (higher growth rate) strains.

Figure 4.6: Stable states found with the temporal method with duration at each oxygen diffusivity value of 1e6 hours. The left column of panels refers to the stable state when oxygen diffusivity is at first high (favouring the oxic state). The right column refers to when oxygen diffusivity is at first low (favouring the anoxic state). Points (which sometimes are so close as to appear as thick lines) show the stable state at the end of each period of constant oxygen diffusivity. Thin lines are for visualisation only, and join the points. Grey arrows show the oxygen diffusivity at which the system flips, and the direction of the flip. The length of the grey arrows is without meaning.

Figure 4.7: Stable states found with the temporal method with duration at each oxygen diffusivity value of 2e6 hours. All other aspects are as in Figure 4.6.

5 SB diversity effects

5.1 SB diversity effect on the anoxic to oxic transition

We showed above that wait times of 1e5 (Figure 4.2) and 1e4 (Figure 4.4) produce qualitatively and often quantitatively very similar results to those for wait times of 1e6 (Figure 4.1). There is, however, a notable difference when there is variation in only SB, or in CB-SB with wait times of 1e4 (Figure 4.4). Here, moderate increases in diversity have the same effect as with wait times of 1e6. However, larger increases in diversity suddenly cause a large reduction in resilience of the anoxic-oxic transition. Figure 5.1 gives an example of the stable states when diversity is moderate – increases in oxygen diffusivity cause strain replacement in the SB to the more tolerant strains, which then delay the shift to the oxic state. Suprisingly, when there is just a little more diversity in the SB, then the anoxic to oxic transition happens earlier, i.e. there is a large decrease in resilience (Figure 5.2).

Figure 5.1: Stable state of the system with moderate level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure 4.6.

Figure 5.2: Stable state of the system with higher than moderate level of diversity in only SB, and no diversity in the other two groups. All other aspects are as in Figure 5.1.

We believe that this phenomenon is caused by an effect that we explain with reference to an analogy of using stepping stones to cross a river. Increasing diversity in our simulations is akin to increasing the width of the river, but keeping the number of stepping stones (strains) constant. So long as the river is not too wide, we can still cross it, but once the width of the river is such that the stepping stones are too far apart, we can’t cross the river. And this happens suddenly, the width increases just a tiny bit more, and we suddenly can’t cross the river when we could before.

In the simulated microbial ecosystem, the strains are the stepping stones, and it is harder (takes longer) for strains to replace each other when there is greater distance (in trait space) between them. When they can no longer replace each other in the time available, the system behaves as if only the least tolerant, highest growth rate strain (lightest shade of colour) is present. In the case illustrated above, this leads to much lower resilience of the anoxic system state.

This phenomenon is dependent on having static trait values for the strains. If strain traits were dynamic (e.g. due to evolution), then we expect this result to not occur, and rather for the effects of diversity with two or three strains to be the same as those with six or nine.

5.2 SB diversity effect on the oxic to anoxic transition

When the environment ameliorates for the sulfate-reducing bacteria (i.e. trajectory of decreasing oxygen diffusivity), there is only limited evidence of strain replacement (in final state) along the oxygen diffusivity gradient – either the system is oxic or the least tolerant strain dominates. And yet the switch to anoxic state occurs at higher oxygen diffusivity than when there is no diversity, which suggests there is some role of the more tolerant strains. We here show that the more tolerant strains do play a role, but only a transient one. The following dynamics are for the system starting oxic, and with a value of oxygen diffusivity (log10(oxygen diffusivity) = -4.8). When there is no diversity within functional groups (i.e. two strains that are identical) the system remains oxic (Figure 5.3). When there is diversity (i.e. the two strains differ) then the system switches to the anoxic state (Figure 5.4). In these examples there is only diversity in the two groups of sulfur bacteria. The most tolerant strain does at first grow, but is then outcompeted by less tolerant strains as the concentration of the inhibiting substrate (oxygen) declines.

Figure 5.3: Dynamics of the system with no diversity in any functional group, and log10(oxygen diffusivity) = -4.8. The system remains oxic when starting in the oxic state. Although two strains are shown in the legend, they have identical parameters, and so their dynamics are identical and only one is visible in the dynamics.

Figure 5.4: Dynamics of the system with diversity in the two groups of sulfur bacteria and otherwise exactly the same conditions as in Figure 5.3. The system flips from oxic to anoxic, whereas it did not when there was no diversity.

6 PB diversity effects

Increases in the diversity within the PB functional group leads to lower resilience of the anoxic state (Figure 4.1, Panel PB, blue line). For example, a small amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately \(log_{10}\)(oxygen diffusivity) = -2.7 (Figure 6.1). Whereas a moderate amount of variation in only the PB functional group, and no variation in other functional groups, results in transition from anoxic to oxic at approximately \(log_{10}\)(oxygen diffusivity) = -3.34 (Figure 6.2).

Figure 6.1: Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a low amount of variation in only the PB functional group. The transition from anoxic to oxic occurs at approximately \(log_{10}\)(oxygen diffusivity) = -2.7. All other aspects are the same as in Figure 3.1.

Figure 6.2: Right-hand panels show the transition from anoxic to oxic state when oxygen diffusivity decreases. The results here are for a moderate amount of variation in only the PB functional group. The transition from anoxic to oxic occurs at approximately \(log_{10}\)(oxygen diffusivity) = -3.34. All other aspects are the same as in Figure 3.1.

This result is different to the effect of diversity within the CB and SB functional groups, which both cause increases in resilience of the state that they dominate (Figure 4.1). Why does diversity in the PB functional group reduce the resilience of the state it dominates?

The counter-intuitive effect of diversity in the phototrophic sulfur bacteria (PB) possibly results from their positive effects on conditions and processes that favour the oxic state. The phototrophic sulfur bacteria oxidise reduced sulfur to oxidised sulfur. Since reduced sulfur inhibits the cyanobacteria (CB), PB can have a positive effect on CB via their consumption of reduced sulfur. Hence, high growth rates of PB will have a positive effect on CB, and therefore favour the oxic state. We would expect to see that higher diversity of PB is associated with lower concentrations of reduced sulfur and higher concentrations of oxidised sulfur before the transition.

Another pathway by which a higher growth rate of PB can favour the oxic state is by indirectly increasing the concentration of oxygen, thus inhibiting themselves and sulfate-reducing bacteria. This effect occurs because PB consume reduced sulfur, such that less reduced sulfur is available for abiotic oxidation. Consequently, a larger amount of oxygen remains in the environment and inhibits sulfur bacteria. We would expect to see that higher diversity of PB is associated with higher concentrations of oxygen before the transition.

Among these two pathways, effects of phototrophic sulfur bacteria via abiotic oxidation of sulfide are probably of subordinate importance in natural environments because abiotic oxidation rates of sulfide are very small compared to biotic oxidation rates (Luther et al 2011).

There is, however, at least one pathway by which PB can favour the anoxic state: PB produce oxidised sulfur which is consumed by SB, and SB produce reduced sulfur, which inhibits CB.

Looking at the strain replacement that occurs across the oxygen diffusivity gradient, we see that as oxygen diffusivity increases, the more oxygen tolerant strains of PB replace the less tolerant ones because they have higher realised growth rates under the prevailing environmental conditions. Greater PB diversity leads to the availability of strains with higher oxygen tolerance and thus even higher realised growth rates, which via the pathways described above, can favour the oxic state.

Substrate concentrations close to the anoxic-oxic transition in simulations with low and high PB diversity, respectively, corroborate the expected effects (Figure 6.3). As the transition is approached, higher PB diversity leads to higher concentration of oxygen, higher concentration of oxidised sulfur (SO), and lower concentration of reduced sulfur (SR), than a simulation with lower PB diversity. All these differences in substrate concentrations are in the direction of favouring the oxic state, resulting in the earlier transition to the oxic state.

Figure 6.3: Concentrations of three substrates (O: oxygen; SO: oxidised sulfur; SR: reduced sulfur) in a low PB diversity simulation (x-axis) versus substrate concentration in a high PB diversity simulation (y-axis). The figures show only a small range of oxygen diffusivity on the trajectory of increasing oxygen diffusivity. The black line is the 1:1 line.

7 Effect of number of strains

In the main text, we present results of simulations with nine strains per functional group. Here we also explore the effect of number of strains (2, 3, 6, 9) on the position of the two tipping points (Figure 7.1). At low to moderate amounts of trait variation, the number of strains did not affect the position of the tipping points. However, at high amounts of trait variation, resilience was lower when the functional groups had only two or three strains rather than six or nine. This effect was particularly pronounced when there was high variation in only SB or in SB-CB: here, simulations with two or three strains led to dramatically lower resilience of the anoxic state than simulations with six or nine strains. We believe that this effect is due to the greater distance in trait space between strains in simulations with two or three strains than in simulations with more strains (also see Section 5).

Figure 7.1: Effects of number of strains and of diversity among strains.

8 Is it really functional diversity?

Reviewer 1: “This made me wonder whether the same (or at least a very similar) results regarding the change in resilience of the system states could also be obtained with just a single strain per functional group, whose tolerance to the inhibiting substance is systematically increased. This could be tested fairly easily by re-running the simulations behind Fig. 3 with just the most tolerant (and slowest growing) strain of each functional group present.”

We made new simulations with communities composed of functional groups all containing nine strains just as before, but then only retained in each functional group the one strain with the greatest tolerance. E.g., in a simulation with variation in the only the SB group (and no variation in the CB or PB groups) then the CB and PB groups (as before) each had only one strain with the mean tolerance value (which is also the maximum there), whereas the SB group contained one strain with higher tolerance. We also performed simulations in which we retained 9 strains in each functional group, but set the tolerance of all nine to the maximum tolerance value. Results were the same as those presented below

For most of the explored combinations and values of functional diversity the results are the same (Figure 8.1). There are differences when there is high functional diversity / high tolerance in the SB, in the CB-SB, and in the CB-SB-PB combinations. (These are also the situations in which we had already found that number of strains can be important (Fig. 7.1 in the original version of the Supplement).) We look at these more closely below.

Figure 8.1: Effects of strain variation on position of tipping points compared between simulations with all (nine) strains present and simulations where only the most tolerant strain is present.

8.1 SB only, high diversity

With only the most tolerance SB strain, the position of the oxic to anoxic tipping point at high SB diversity appears to be at very low oxygen diffusivity, compared to when all nine strains are present (Figure 8.1). A corresponding simulations from which the position of the tipping points are calculated are shown in Figures 8.2 (nine strains) and 8.3 (most tolerant strain only). We see that the effect of oxygen diffusivity on the stable state is quite different in the case of only the most tolerant strain only–at low oxygen diffusivity there is variation around a stable state with all three functional groups coexisting–this never occured when there were multiple strains varying in tolerance, or when there was only one strain with mean tolerance. Furthermore, the state of the system is independent of the direction of change in oxygen diffusivity. I.e. there is no region of bistability (no alternate stable states).

The algorithm used to detect state changes is working incorrectly in this specific situation. Here is the description of the algorithm from the main text: “We calculated the tipping point position of the oxic to anoxic (anoxic-oxic transition as the minimum (maximum) level of oxygen diffusivity at which the total density of sulfate-reducing bacteria was greater than 1000 abundance units different between the increasing and decreasing oxygen diffusivity simulation.” This is not of great concern, however, as there are no state transitions to detect.

We now wish to return to an aspect of the reviewer comment: It is correct that in Fig. 4 when oxygen diffusivity increased, the two dominant functional groups (the sulfur bacteria groups) are dominated by the most tolerant strain before the state transition. And when oxygen diffusivity decreased, the dominant functional group (the cyanobacteria) is dominated by the most tolerant strain before the state transition. It is worth noting, however, that the transitions that occur are to states dominated by the least tolerant-highest growth rate strain. This raises the possibility that less tolerant strains may influence the results.

Indeed, in the SB case, absence of the high growth rate SB strain precludes the formation of a stable state that does not include cyanobacteria. It seems that the most tolerant SB strain does not have a high enough growth rate to itself exclude the cyanobacteria. Cyanobacteria are therefore always present, regardless of the oxygen diffusivity. (The stable states of the system then display a complex relationship with oxygen diffusivity that our algorithm for detecting tipping points is not well suited to, hence the rather odd positions of the tipping points.)

Figure 8.2: Stable states found with nine strains of SB and high diversity. All other aspects are the same as in Figure 3.1.

Figure 8.3: Stable states found with the single most tolerant SB and high diversity. All other aspects are the same as in Figure 3.1.

8.1.1 What causes variation in response to oxygen diffusivity (SB)?

In the ecosystem shown in Figure 8.3, where there was quite high tolerance SB strain, there was large variation in state variables along the oxygen diffusivity gradient at values less than about -4.5. To examine what might be causing this, we simulated the ecosystem at a constant oxygen diffusivity of \(log_{10}\) oxygen diffusivity of -6 (Figure 8.4). We see that the dynamics are qualitatively different from other cases in which there is a stable state reached in such a simulation. It appears that the high tolerance SB strain is able to initiate a transition to an anoxic state, but is not able to create a stable anoxic stable state. This dynamical situation was not presented in the ecosystem in which high tolerance and low tolerance strains were present.

## [1] 2e+05

Figure 8.4: Temporal dynamics of an ecosystem with a single SB strain with high tolerance, the same one as in Figure 8.3 simulated with constant $log_{10}(oxygen diffusivity) of -6.

8.2 CB-SB only, high diversity

As in the SB only case, the absence of lower tolerance-higher growth rate strains changes the system to a parameter space where it does not display alternate states, at least within the range of oxygen diffusivities simulated.

Figure 8.5: Stable states found with nine strains of CB and SB and high diversity. All other aspects are the same as in Figure 3.1.

Figure 8.6: Stable states found with the single most tolerant CB strain and SB strain and high diversity. All other aspects are the same as in Figure 3.1.

8.2.1 What causes variation in response to oxygen diffusivity (CB-SB)?

In the ecosystem shown in Figure 8.6, where there was quite high tolerance CB strain and SB strain, there was large variation in state variables along the oxygen diffusivity gradient at values less than about -6.5. To examine what might be causing this, we simulated the ecosystem at a constant oxygen diffusivity of \(log_{10}\) oxygen diffusivity of -6 (Figure 8.7). We see that the dynamics are qualitatively different from other cases in which there is a stable state reached in such a simulation. It appears that the high tolerance SB strain is able to initiate a transition to an anoxic state, but is not able to create a stable anoxic stable state.

## [1] 2e+05

Figure 8.7: Temporal dynamics of an ecosystem with a single CB strain with high tolerance and a single SB strain with high tolerance, the same one as in Figure 8.6 simulated with constant $log_{10}(oxygen diffusivity) of -6.

8.3 CB-SB-PB, high diversity

In the CB-SB-PB case when only the most tolerant strain of all functional groups is present, the resilience of the oxic state is greater than when variation in strain tolerance is present. Looked at from another perspective, comparing CB-SB and CB-SB-PB both with only the most tolerant strains, we see that presence of a high tolerance PB strain (in CB-SB-PB) allows the formation of an anoxic cyanobacteria-free stable state (whereas there was no such state in CB-SB). We interpret this as presence of a high tolerance PB strain stabilising the anoxic state, perhaps by increasing the availability of oxidised sulfur which is a resource for SB.

Figure 8.8: Stable states found with nine strains of CB-SB-PB and high diversity. All other aspects are the same as in Figure 3.1.

Figure 8.9: Stable states found with the single most tolerant CB strain, SB strain, and PB strain, and high diversity. All other aspects are the same as in Figure 3.1.

8.3.1 Why only at high diversity?

Why do these effects only occur at high diversities? The likely reason is that at high diversities then the most tolerant strain has particularly low growth rate, and so it may have the tolerance required to persist, it does not have the growth rate required to dominate.

9 Finding the separatrix

In the original publication (Bush et al 2017) Figure 3a shows a separatrix in the cyanobacteria abundance - oxygen diffusivity dimensions of the system. If the system has initial cyanobacterial abundance greater than the value at the separatrix then the cyanobacterial abundance increases, and the system goes to the oxic stable state. Conversely, if the system has initial cyanobacterial abundance lower than the value at the separatrix then the cyanobacterial abundance decreases, and the system goes to the anoxic stable state.

It is critical to note, however, that the position of the separatrix depends on both the parameters of the system, and the initial conditions of variables other than cyanobacterial abundance. If the system were started with higher initial concentration, then the position of the separatrix would change (as would the values of oxygen diffusivity at which the system tipped from anoxic to oxic, and from oxic to anoxic).

This issue of dependence of the position of the separatrix on the value of initial condition of multiple state variables has particular significance when there are multiple strains per functional group. Then one cannot just vary initial total cyanobacterial abundance (as was done to create the separatrix on Figure 3a in the original publication) but one must also define how the initial abundance of each of the strains varies. One might choose to make each strain have equal initial abundance, and to vary the total (as we did in our implementation of the replication method in the main article).

We could, however, and arguably should, give the most tolerant strain of cyanobacteria higher initial relative abundance when total initial cyanobacterial abundance was low, and the least tolerant strain higher initial relative abundance when total cyanobacterial abundance was high. And we would need to specify the total and relative abundance of all cyanobacteria (and sulfur bactria) strains. Thus, with a system as complex as the 9 strain system, which has 31 state variables the separatrix is very high dimensional, and evaluating (and plotting) it in one dimension requires many assumptions. We chose to not attempt to find a meaningful set of assumptions in this study, though not because we think it would be unprofitable to do so, but rather because we did not want to risk distracting attention from findings we consider to be significant enough in their own right.

10 (Non-)additivity

We investigated if the effects of trait variation of individual functional groups were additive. We compared the tipping point positions as determined by our simulations with the expected position of the respective tipping point if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes (i.e. the difference in tipping point position between the no-diversity scenario and a given amount of trait variation) on a linear scale.

Figure 10.1: Assessing the additivity or non-additivity of the effects of trait variation. Solid black lines depict the position of the two tipping points at different levels of trait variation, as simulated with the model. The colored, dashed lines show the expected position of tipping points if effects of diversity in different functional groups were additive. The expected tipping point was calculated by adding the effect sizes of trait variation of different functional groups on a linear scale. For example, in (d) the effect size of variation in only CB was added to the effect size of variation in only SB. In (g), the effect size of variation in only CB was added to the effect size of simultaneous variation in both SB and PB. (CB: cyanobacteria, SB: sulfate-reducing bacteria, PB: phototrophic sulfur bacteria). The dotted vertical lines at -8 and 0 on the x-axis show the extent of the oxygen diffusivity gradient investigated in the simulations.

When two functional groups varied in traits, the effects of trait variation were additive when one of the groups had no diversity effect on a tipping point but non-additive when both groups had diversity effects on a tipping point (Fig. 10.1 a-f). High trait variation in cyanobacteria and sulfate-reducing bacteria caused the oxic-anoxic tipping point to occur at lower oxygen diffusivity than expected assuming additivity of the individual effects of diversity in these two groups (Fig. 10.1d). When the two groups of sulfur bacteria varied in traits, the anoxic-oxic tipping point occurred at lower oxygen diffusivity than expected by additivity (Fig. 10.1f).

When all three groups varied in traits, only the diversity effect in phototrophic sulfur bacteria on the oxic-anoxic tipping point was additive; all other effects were non-additive, especially at higher trait variation (Fig. 10.1g-i). Both tipping points occurred at lower levels of oxygen diffusivity than expected by additivity. Interestingly, diversity in the cyanobacteria moved the shift to the oxic state to lower levels of oxygen diffusivity although independent variation in the cyanobacteria had no effect on this tipping point.

11 References

Bush, T., Diao, M., Allen, R.J., Sinnige, R., Muyzer, G. & Huisman, J. (2017). Oxic-anoxic regime shifts mediated by feedbacks between biogeochemical processes and microbial community dynamics. Nat. Commun., 8, 789.

Luther, G.W., Findlay, A., MacDonald, D., Owings, S., Hanson, T., Beinart, R., et al. (2011). Thermodynamics and kinetics of sulfide oxidation by oxygen: a look at inorganically controlled reactions and biologically mediated processes in the environment. Front. Microbiol., 2, 62.